---
title: "Code"
format: html
---

```{r}
if (!require("pacman")) install.packages("pacman")
pacman::p_load(readr,dplyr,ggplot2, sandwich, forcats, lmtest, stringr, broom, flextable, usdata)
```



```{r}
data <- read_csv("tweet_level_data.csv")
```


# Comparison between mandate-related and non-mandate-related tweets
## Freedom-related words

Ha: mandate-related tweets were more likely to contain freedom related words 

proportional test
```{r}
# z value
(mean(data$pandemic_freedom[data$pred==1])-mean(data$pandemic_freedom[data$pred==0]))/sqrt(var(data$pandemic_freedom[data$pred==1])/sum(data$pred==1)+var(data$pandemic_freedom[data$pred==0])/sum(data$pred==0))
# p-value
pnorm((mean(data$pandemic_freedom[data$pred==1])-mean(data$pandemic_freedom[data$pred==0]))/sqrt(var(data$pandemic_freedom[data$pred==1])/sum(data$pred==1)+var(data$pandemic_freedom[data$pred==0])/sum(data$pred==0)),lower.tail = FALSE)
```

## Negativity

Ha: Mandate-related tweets were more negative in VADER negative scores

Difference in means by the central limit theorem
```{r}
# z value
(mean(data$vader_neg[data$pred==1])-mean(data$vader_neg[data$pred==0]))/sqrt(var(data$vader_neg[data$pred==1])/sum(data$pred==1)+var(data$vader_neg[data$pred==0])/sum(data$pred==0))
# p-value
pnorm((mean(data$vader_neg[data$pred==1])-mean(data$vader_neg[data$pred==0]))/sqrt(var(data$vader_neg[data$pred==1])/sum(data$pred==1)+var(data$vader_neg[data$pred==0])/sum(data$pred==0)),lower.tail = FALSE)
```

Wilcoxon signed-rank test
```{r}
wilcox.test(data$vader_neg[data$pred==1], data$vader_neg[data$pred==0], alternative = "greater")
```

Ha: Mandate-related tweets were more likely to be negative (classified by TweetNLP)

```{r}
# z value
(mean(data$nlp_neg_label[data$pred==1])-mean(data$nlp_neg_label[data$pred==0]))/sqrt(var(data$nlp_neg_label[data$pred==1])/sum(data$pred==1)+var(data$nlp_neg_label[data$pred==0])/sum(data$pred==0))
# p-value
pnorm((mean(data$nlp_neg_label[data$pred==1])-mean(data$nlp_neg_label[data$pred==0]))/sqrt(var(data$nlp_neg_label[data$pred==1])/sum(data$pred==1)+var(data$nlp_neg_label[data$pred==0])/sum(data$pred==0)),lower.tail = FALSE)
```

## Anger

Ha: Mandate-related tweets were angrier (measured by NRC Lexicon)

Difference in means by the central limit theorem
```{r}
(mean(data$NRCLex_anger[data$pred==1])-mean(data$NRCLex_anger[data$pred==0]))/sqrt(var(data$NRCLex_anger[data$pred==1])/sum(data$pred==1)+var(data$NRCLex_anger[data$pred==0])/sum(data$pred==0))
pnorm((mean(data$NRCLex_anger[data$pred==1])-mean(data$NRCLex_anger[data$pred==0]))/sqrt(var(data$NRCLex_anger[data$pred==1])/sum(data$pred==1)+var(data$NRCLex_anger[data$pred==0])/sum(data$pred==0)),lower.tail = FALSE)
```

Wilcoxon signed-rank test
```{r}
wilcox.test(data$NRCLex_anger[data$pred==1], data$NRCLex_anger[data$pred==0], alternative = "greater")
```

Ha: Mandate-related tweets were more likely to be angry (classified by TweetNLP)

```{r}
# z value
(mean(data$nlp_anger[data$pred==1])-mean(data$nlp_anger[data$pred==0]))/sqrt(var(data$nlp_anger[data$pred==1])/sum(data$pred==1)+var(data$nlp_anger[data$pred==0])/sum(data$pred==0))
# p-value
pnorm((mean(data$nlp_anger[data$pred==1])-mean(data$nlp_anger[data$pred==0]))/sqrt(var(data$nlp_anger[data$pred==1])/sum(data$pred==1)+var(data$nlp_anger[data$pred==0])/sum(data$pred==0)),lower.tail = FALSE)
```

## Figure 2
```{r}
all_avg_report <- rbind(data %>% group_by(pred)%>%
  summarise(y=mean(pandemic_freedom*100,na.rm=TRUE),
            se=sd(pandemic_freedom*100,na.rm=TRUE)/sqrt(n()))%>%
  mutate(var="Freedom"),
data %>% group_by(pred)%>%
  summarise(y=mean(vader_neg*100,na.rm=TRUE),
            se=sd(vader_neg*100,na.rm=TRUE)/sqrt(n()))%>%
  mutate(var="Rule-Based\nNegativity"),
data %>% group_by(pred)%>%
  summarise(y=mean(nlp_neg_label*100,na.rm=TRUE),
            se=sd(nlp_neg_label*100,na.rm=TRUE)/sqrt(n()))%>%
  mutate(var="ML Negativity"),
data %>% group_by(pred)%>%
  summarise(y=mean(NRCLex_anger*100,na.rm=TRUE),
            se=sd(NRCLex_anger*100,na.rm=TRUE)/sqrt(n()))%>%
  mutate(var="Dictionary-Based\nAnger"),
data %>% group_by(pred)%>%
  summarise(y=mean(nlp_anger*100,na.rm=TRUE),
            se=sd(nlp_anger*100,na.rm=TRUE)/sqrt(n()))%>%
  mutate(var="ML Anger"))%>%
  mutate(var=factor(var,levels=c("Freedom","Rule-Based\nNegativity","ML Negativity", "Dictionary-Based\nAnger", "ML Anger")))
```

```{r}
ggplot(all_avg_report,aes(x=factor(pred,label=c("No","Yes")), y=y))+
  geom_bar(stat="identity",position = position_dodge())+
  geom_errorbar(aes(ymin=y+qnorm(0.025)*se, ymax=y+qnorm(0.975)*se), width=.2,position=position_dodge(.9)) +
  geom_text(aes(label = round(y,3)), vjust = 1.5, colour = "white",size = 3)+
  scale_x_discrete(name="Mandate-related")+
  scale_y_continuous(name="Values (0-100)")+
  facet_wrap(~var,scales = "free_y",nrow=2)+
  theme_bw(base_size = 10)+
  theme(legend.position="bottom")
ggsave("p_results_by_mandate.tiff",device="tiff",width=1280,height=1280,units = "px")
```

## Figure S1
```{r}
all_avg_sentiment <- rbind(
data %>% group_by(pred)%>%
  summarise(y=mean(vader_pos*100,na.rm=TRUE),
            se=sd(vader_pos*100,na.rm=TRUE)/sqrt(n()))%>%
  mutate(var="Positive"),
data %>% group_by(pred)%>%
  summarise(y=mean(vader_compound*100,na.rm=TRUE),
            se=sd(vader_compound*100,na.rm=TRUE)/sqrt(n()))%>%
  mutate(var="Compound"),
data %>% group_by(pred)%>%
  summarise(y=mean(vader_neu*100,na.rm=TRUE),
            se=sd(vader_neu*100,na.rm=TRUE)/sqrt(n()))%>%
  mutate(var="Neutral")
)%>%
  mutate(var=factor(var,levels=c("Positive","Neutral","Compound")))
```

```{r}
ggplot(all_avg_sentiment,aes(x=factor(pred,label=c("No","Yes")), y=y))+
  geom_bar(stat="identity",position = position_dodge())+
  geom_errorbar(aes(ymin=y+qnorm(0.025)*se, ymax=y+qnorm(0.975)*se), width=.2,position=position_dodge(.9)) +
  geom_text(aes(label = round(y,3)), vjust = 1.5, colour = "grey",size = 4)+
  scale_x_discrete(name="Mandate-related")+
  scale_y_continuous(name="Values")+
  facet_wrap(~var,scales = "free_y",nrow=1)+
  theme_bw(base_size = 15)+
  theme(legend.position="bottom")
```

## Figure S2
```{r}
all_avg_popularity <- rbind(data %>% group_by(pred)%>%
  summarise(y=mean(retweet_count,na.rm=TRUE),
            se=sd(retweet_count,na.rm=TRUE)/sqrt(n()))%>%
  mutate(var="Retweet"),
data %>% group_by(pred)%>%
  summarise(y=mean(reply_count,na.rm=TRUE),
            se=sd(reply_count,na.rm=TRUE)/sqrt(n()))%>%
  mutate(var="Reply"),
data %>% group_by(pred)%>%
  summarise(y=mean(like_count,na.rm=TRUE),
            se=sd(like_count,na.rm=TRUE)/sqrt(n()))%>%
  mutate(var="Like"))
```

```{r}
ggplot(all_avg_popularity,aes(x=factor(pred,label=c("No","Yes")), y=y))+
  geom_bar(stat="identity",position = position_dodge())+
  geom_errorbar(aes(ymin=y+qnorm(0.025)*se, ymax=y+qnorm(0.975)*se), width=.2,position=position_dodge(.9)) +
  geom_text(aes(label = round(y,3)), vjust = 1.5, colour = "white",size = 4)+
  scale_x_discrete(name="Mandate-related")+
  scale_y_continuous(name="Counts")+
  facet_wrap(~var,scales = "free_y")+
  theme_bw(base_size = 15)+
  theme(legend.position="bottom")
```

# Figure 1

```{r}
data_ts2 <- data %>%
  filter(state %in% state_stats$state)%>%
  group_by(date) %>%
  summarise(p_man=sum(pred==1)/n(),
            nlp_anger=mean(nlp_anger,na.rm = TRUE))
  
```

```{r}
data_ts2 <- data %>%
  filter(state %in% state_stats$state)%>%
  group_by(date) %>%
  summarise(p_man=sum(pred==1)/n(),
            nlp_anger=mean(nlp_anger,na.rm = TRUE))
  
```


```{r}
ggplot(data_ts2) + 
  geom_vline(aes(xintercept=date),
             data_ts2%>%filter(date%in%c(as.Date("2021-07-26"),
                                        as.Date("2021-08-03"),
                                        as.Date("2021-09-09"),
                                        as.Date("2021-10-12"),
                                        as.Date("2021-11-04"),
                                        as.Date("2022-01-13"))),colour="black",linetype="dashed",linewidth=0.2)+
  geom_line(aes(x=date,y=p_man*100,color="Mandate",linetype="Mandate"),linewidth=0.4)+
  geom_line(aes(x=date,y=nlp_anger*100,color="Anger",linetype="Anger"),linewidth=0.4)+
  geom_label(aes(x=as.Date("2021-07-26"),y=3,label="Jul 26, 2021\n California statewide\n vaccine mandate"),size=2,nudge_x = -3)+
  geom_label(aes(x=as.Date("2021-08-03"),y=13,label="Aug 3, 2021\n NYC indoor\n vaccination requirements"),size=2)+
  geom_label(aes(x=as.Date("2021-09-09"),y=3,label="Sep 9, 2021\n White House\n announcement"),size=2)+
  geom_label(aes(x=as.Date("2021-10-12"),y=15,label="Oct 12, 2021\n Proposed federal\n mandate released"),size=2)+
  geom_label(aes(x=as.Date("2021-11-04"),y=5,label="Nov 4, 2021\n Federal mandates\n issued"),size=2)+
  geom_label(aes(x=as.Date("2022-01-13"),y=5,label="Jan 13, 2022\n Supreme court\n rulings"),size=2)+
  scale_x_date(date_breaks="1 month",date_labels = "%b%Y")+
  scale_y_continuous(name="Percentage of Tweets",limits=c(0,53))+
  scale_color_manual(name="",breaks=c("Mandate","Anger"),values=c("Mandate"="#547c99","Anger"="#93003a"))+
  scale_linetype_manual(name="",breaks=c("Mandate","Anger"),values=c("Mandate"="solid","Anger"="longdash"))+
  theme_bw(base_size = 8)+
  theme(legend.position="bottom")
ggsave("s_p_man_neg_ts.tiff",device="tiff",width=1380,height=850,units = "px")

```

# State-date panel data analysis

```{r}
data_panel <- read_csv("panel_data.csv")
```


DV: the percentage of vaccine-related tweets containing freedom-related words
```{r}
lm_word_freedom_twfe1 <- lm(word_freedom*100~-1+I(p_man*10)+log(n_100k.vac)+
                               log_new_case_7ma+log_new_death_7ma+log_new_dose_7ma+
              I(Series_Complete_Pop_Pct/10)+I(Series_Complete_Pop_Pct_sd/10)+
                  as.factor(date)+as.factor(state),
                       data=data_panel)
vcov_word_freedom_twfe1 <- vcovPL(lm_word_freedom_twfe1, 
                                          cluster = data_panel$state,
                                          order.by = data_panel$date)
output_lm_word_freedom_twfe1 <- coeftest(lm_word_freedom_twfe1, vcov = vcov_word_freedom_twfe1)
output_lm_word_freedom_twfe1<-output_lm_word_freedom_twfe1[str_detect(row.names(output_lm_word_freedom_twfe1),"^(?!as)"),]
output_lm_word_freedom_twfe1
```

Figure 3
```{r}
data.frame(dv="Freedom words%",
           iv=factor(names(output_lm_word_freedom_twfe1[,1]),
                         levels = names(output_lm_word_freedom_twfe1[,1]),
                         labels = c("Mandate tweets%\n(Every 10%)",
                                    "log n\nvaccine tweets",
                                    "log new case",
                                    "log new death",
                                    "log new dose",
                                    "Vaccine coverage",
                                    "Vaccine coverage sd")),
           estimate=output_lm_word_freedom_twfe1[,1],
           lb95=output_lm_word_freedom_twfe1[,1]+qnorm(0.025)*output_lm_word_freedom_twfe1[,2],
           ub95=output_lm_word_freedom_twfe1[,1]+qnorm(0.975)*output_lm_word_freedom_twfe1[,2],
           lb90=output_lm_word_freedom_twfe1[,1]+qnorm(0.05)*output_lm_word_freedom_twfe1[,2],
           ub90=output_lm_word_freedom_twfe1[,1]+qnorm(0.95)*output_lm_word_freedom_twfe1[,2])%>%
  mutate(iv=fct_rev(iv))%>%
  ggplot(aes(iv, estimate))+
  geom_point()+
  geom_text(aes(label=round(estimate,3)),nudge_x = 0.4,size=2.5)+
  geom_errorbar(aes(ymin = lb95, ymax = ub95), width=.2)+
  geom_errorbar(aes(ymin = lb90, ymax = ub90), width=.4)+
  geom_hline(yintercept = 0,
             colour = gray(1/2), lty = 2)+
  scale_x_discrete(name="Independent Variables")+
  scale_y_continuous(name="Effects")+
  coord_flip()+
  theme_bw(base_size=12)  
ggsave("p_panel_freedom.tiff",width=1280,height=720,units = "px")
```

DV: average VADER negative sentiment score 
```{r}
lm_vader_neg.vac_twfe1 <- lm(vader_neg*100~-1+I(p_man*10)+log(n_100k.vac)+
                               log_new_case_7ma+
                               log_new_death_7ma+
                               log_new_dose_7ma+
              I(Series_Complete_Pop_Pct/10)+I(Series_Complete_Pop_Pct_sd/10)+
                  as.factor(date)+as.factor(state),
                       data=data_panel)
vcov_vader_neg.vac_twfe1 <- vcovPL(lm_vader_neg.vac_twfe1, 
                                          cluster = data_panel$state,
                                          order.by = data_panel$date)
output_lm_vader_neg.vac_twfe1 <- coeftest(lm_vader_neg.vac_twfe1,
                                   vcov = vcov_vader_neg.vac_twfe1)

output_lm_vader_neg.vac_twfe1<-output_lm_vader_neg.vac_twfe1[str_detect(row.names(output_lm_vader_neg.vac_twfe1),"^(?!as)"),]
output_lm_vader_neg.vac_twfe1

```

DV: percentage of negative tweets classified by TweetNLP 
```{r}
lm_nlp_neg.vac_twfe1 <- lm(nlp_neg_label*100~-1+I(p_man*10)+log(n_100k.vac)+
                               log_new_case_7ma+
                               log_new_death_7ma+
                               log_new_dose_7ma+
              I(Series_Complete_Pop_Pct/10)+I(Series_Complete_Pop_Pct_sd/10)+
                  as.factor(date)+as.factor(state),
                       data=data_panel)
vcov_nlp_neg.vac_twfe1 <- vcovPL(lm_nlp_neg.vac_twfe1, 
                                          cluster = data_panel$state,
                                          order.by = data_panel$date)
output_lm_nlp_neg.vac_twfe1 <- coeftest(lm_nlp_neg.vac_twfe1,
                                   vcov = vcov_nlp_neg.vac_twfe1)

output_lm_vader_neg.vac_twfe1<-output_lm_nlp_neg.vac_twfe1[str_detect(row.names(output_lm_nlp_neg.vac_twfe1),"^(?!as)"),]
output_lm_vader_neg.vac_twfe1

```

DV: average percentage of anger words from the NRC Lexicon
```{r}
lm_NRC_anger.vac_twfe1 <- lm(NRC_anger*100~-1+I(p_man*10)+log(n_100k.vac)+
                               log_new_case_7ma+log_new_death_7ma+log_new_dose_7ma+
              I(Series_Complete_Pop_Pct/10)+I(Series_Complete_Pop_Pct_sd/10)+
                  as.factor(date)+as.factor(state),
                       data=data_panel)
vcov_NRC_anger.vac_twfe1 <- vcovPL(lm_NRC_anger.vac_twfe1, 
                                          cluster = data_panel$state,
                                          order.by = data_panel$date)

output_lm_NRC_anger.vac_twfe1 <- coeftest(lm_NRC_anger.vac_twfe1, vcov = vcov_NRC_anger.vac_twfe1)
output_lm_NRC_anger.vac_twfe1<-output_lm_NRC_anger.vac_twfe1[str_detect(row.names(output_lm_NRC_anger.vac_twfe1),"^(?!as)"),]
output_lm_NRC_anger.vac_twfe1
```

DV: percentage of angry tweets classified by TweetNLP 
```{r}
lm_nlp_anger.vac_twfe1 <- lm(nlp_anger*100~-1+I(p_man*10)+log(n_100k.vac)+
                               log_new_case_7ma+log_new_death_7ma+log_new_dose_7ma+
              I(Series_Complete_Pop_Pct/10)+I(Series_Complete_Pop_Pct_sd/10)+
                  as.factor(date)+as.factor(state),
                       data=data_panel)
vcov_nlp_anger.vac_twfe1 <- vcovPL(lm_nlp_anger.vac_twfe1, 
                                          cluster = data_panel$state,
                                          order.by = data_panel$date)

output_lm_nlp_anger.vac_twfe1 <- coeftest(lm_nlp_anger.vac_twfe1, vcov = vcov_nlp_anger.vac_twfe1)
output_lm_nlp_anger.vac_twfe1<-output_lm_nlp_anger.vac_twfe1[str_detect(row.names(output_lm_nlp_anger.vac_twfe1),"^(?!as)"),]
output_lm_nlp_anger.vac_twfe1
```

Figure 4
```{r}
data.frame(dv="Rule-Based\nNegativity",
           iv=factor(names(output_lm_vader_neg.vac_twfe1[,1]),
                         levels = names(output_lm_vader_neg.vac_twfe1[,1]),
                         labels = c("Mandate tweets%\n(Every 10%)",
                                    "log n\nvaccine tweets",
                                    "log new case",
                                    "log new death",
                                    "log new dose",
                                    "Vaccine coverage",
                                    "Vaccine coverage sd")),
           estimate=output_lm_vader_neg.vac_twfe1[,1],
           lb95=output_lm_vader_neg.vac_twfe1[,1]+qnorm(0.025)*output_lm_vader_neg.vac_twfe1[,2],
           ub95=output_lm_vader_neg.vac_twfe1[,1]+qnorm(0.975)*output_lm_vader_neg.vac_twfe1[,2],
           lb90=output_lm_vader_neg.vac_twfe1[,1]+qnorm(0.05)*output_lm_vader_neg.vac_twfe1[,2],
           ub90=output_lm_vader_neg.vac_twfe1[,1]+qnorm(0.95)*output_lm_vader_neg.vac_twfe1[,2])%>% 
  bind_rows(data.frame(dv="ML-based\nNegativity",
           iv=factor(names(output_lm_nlp_neg.vac_twfe1[,1]),
                         levels = names(output_lm_nlp_neg.vac_twfe1[,1]),
                         labels = c("Mandate tweets%\n(Every 10%)",
                                    "log n\nvaccine tweets",
                                    "log new case",
                                    "log new death",
                                    "log new dose",
                                    "Vaccine coverage",
                                    "Vaccine coverage sd")),
           estimate=output_lm_nlp_neg.vac_twfe1[,1],
           lb95=output_lm_nlp_neg.vac_twfe1[,1]+qnorm(0.025)*output_lm_nlp_neg.vac_twfe1[,2],
           ub95=output_lm_nlp_neg.vac_twfe1[,1]+qnorm(0.975)*output_lm_nlp_neg.vac_twfe1[,2],
           lb90=output_lm_nlp_neg.vac_twfe1[,1]+qnorm(0.05)*output_lm_nlp_neg.vac_twfe1[,2],
           ub90=output_lm_nlp_neg.vac_twfe1[,1]+qnorm(0.95)*output_lm_nlp_neg.vac_twfe1[,2]))%>%
  bind_rows(data.frame(dv="Dictionary-Based\nAnger",
           iv=factor(names(output_lm_NRC_anger.vac_twfe1[,1]),
                         levels = names(output_lm_NRC_anger.vac_twfe1[,1]),
                         labels = c("Mandate tweets%\n(Every 10%)",
                                    "log n\nvaccine tweets",
                                    "log new case",
                                    "log new death",
                                    "log new dose",
                                    "Vaccine coverage",
                                    "Vaccine coverage sd")),
           estimate=output_lm_NRC_anger.vac_twfe1[,1],
           lb95=output_lm_NRC_anger.vac_twfe1[,1]+qnorm(0.025)*output_lm_NRC_anger.vac_twfe1[,2],
           ub95=output_lm_NRC_anger.vac_twfe1[,1]+qnorm(0.975)*output_lm_NRC_anger.vac_twfe1[,2],
           lb90=output_lm_NRC_anger.vac_twfe1[,1]+qnorm(0.05)*output_lm_NRC_anger.vac_twfe1[,2],
           ub90=output_lm_NRC_anger.vac_twfe1[,1]+qnorm(0.95)*output_lm_NRC_anger.vac_twfe1[,2]))%>%
  bind_rows(data.frame(dv="ML-based\nAnger",
           iv=factor(names(output_lm_nlp_anger.vac_twfe1[,1]),
                         levels = names(output_lm_nlp_anger.vac_twfe1[,1]),
                         labels = c("Mandate tweets%\n(Every 10%)",
                                    "log n\nvaccine tweets",
                                    "log new case",
                                    "log new death",
                                    "log new dose",
                                    "Vaccine coverage",
                                    "Vaccine coverage sd")),
           estimate=output_lm_nlp_anger.vac_twfe1[,1],
           lb95=output_lm_nlp_anger.vac_twfe1[,1]+qnorm(0.025)*output_lm_nlp_anger.vac_twfe1[,2],
           ub95=output_lm_nlp_anger.vac_twfe1[,1]+qnorm(0.975)*output_lm_nlp_anger.vac_twfe1[,2],
           lb90=output_lm_nlp_anger.vac_twfe1[,1]+qnorm(0.05)*output_lm_nlp_anger.vac_twfe1[,2],
           ub90=output_lm_nlp_anger.vac_twfe1[,1]+qnorm(0.95)*output_lm_nlp_anger.vac_twfe1[,2]))%>%
  mutate(dv=factor(dv,levels=c("Rule-Based\nNegativity","ML-based\nNegativity","Dictionary-Based\nAnger","ML-based\nAnger")),
         iv=fct_rev(iv))%>%
  ggplot(aes(iv, estimate))+
  geom_point()+
  geom_text(aes(label=round(estimate,3)),nudge_x = 0.4,size=2.5)+
  geom_errorbar(aes(ymin = lb95, ymax = ub95), width=.2, linewidth=0.3)+
  geom_errorbar(aes(ymin = lb90, ymax = ub90), width=.4, linewidth=0.3)+
  geom_hline(yintercept = 0,
             colour = gray(1/2), lty = 2)+
  scale_x_discrete(name="Independent Variables")+
  facet_wrap(~dv,scales="free")+
  scale_y_continuous(name="Effects")+
  coord_flip()+
  theme_bw(base_size=9)
ggsave("p_panel_vax.tiff",width=1920,height=1280,units = "px")
```

DV: average VADER negative sentiment score in tweets about public health officials
```{r}
lm_vader_neg.cdc_twfe1 <- lm(vader_neg.cdc*100~-1+I(p_man*10)+log(n_100k.vac)+
                                 log_new_case_7ma+log_new_death_7ma+log_new_dose_7ma+
              I(Series_Complete_Pop_Pct/10)+I(Series_Complete_Pop_Pct_sd/10)+
                  as.factor(date)+as.factor(state),
                       data=data_panel)
vcov_vader_neg.cdc_twfe1 <- vcovPL(lm_vader_neg.cdc_twfe1, 
                                          cluster = data_panel$state,
                                          order.by = data_panel$date)
output_lm_vader_neg.cdc_twfe1<-coeftest(lm_vader_neg.cdc_twfe1, vcov = vcov_vader_neg.cdc_twfe1)
output_lm_vader_neg.cdc_twfe1<-output_lm_vader_neg.cdc_twfe1[str_detect(row.names(output_lm_vader_neg.cdc_twfe1),"^(?!as)"),]
output_lm_vader_neg.cdc_twfe1
```

DV: percentage of negative tweets classified by TweetNLP in tweets about public health officials
```{r}
lm_nlp_neg.cdc_twfe1 <- lm(nlp_neg_label.cdc*100~-1+I(p_man*10)+log(n_100k.vac)+
                                 log_new_case_7ma+log_new_death_7ma+log_new_dose_7ma+
              I(Series_Complete_Pop_Pct/10)+I(Series_Complete_Pop_Pct_sd/10)+
                  as.factor(date)+as.factor(state),
                       data=data_panel)
vcov_nlp_neg.cdc_twfe1 <- vcovPL(lm_nlp_neg.cdc_twfe1, 
                                          cluster = data_panel$state,
                                          order.by = data_panel$date)
output_lm_nlp_neg.cdc_twfe1<-coeftest(lm_nlp_neg.cdc_twfe1, vcov = vcov_nlp_neg.cdc_twfe1)
output_lm_nlp_neg.cdc_twfe1<-output_lm_nlp_neg.cdc_twfe1[str_detect(row.names(output_lm_nlp_neg.cdc_twfe1),"^(?!as)"),]
output_lm_nlp_neg.cdc_twfe1
```

DV: average percentage of anger words from the NRC Lexicon in tweets about public health officials
```{r}
lm_NRC_anger.cdc_twfe1 <- lm(NRC_anger.cdc*100~-1+I(p_man*10)+log(n_100k.vac)+
                                 log_new_case_7ma+log_new_death_7ma+log_new_dose_7ma+
              I(Series_Complete_Pop_Pct/10)+I(Series_Complete_Pop_Pct_sd/10)+
                  as.factor(date)+as.factor(state),
                       data=data_panel)
vcov_NRC_anger.cdc_twfe1 <- vcovPL(lm_NRC_anger.cdc_twfe1, 
                                          cluster = data_panel$state,
                                          order.by = data_panel$date)

output_lm_NRC_anger.cdc_twfe1<-coeftest(lm_NRC_anger.cdc_twfe1, vcov = vcov_NRC_anger.cdc_twfe1)
output_lm_NRC_anger.cdc_twfe1<-output_lm_NRC_anger.cdc_twfe1[str_detect(row.names(output_lm_NRC_anger.cdc_twfe1),"^(?!as)"),]
output_lm_NRC_anger.cdc_twfe1
```

DV: percentage of angry tweets classified by TweetNLP in tweets about public health officials
```{r}
lm_nlp_anger.cdc_twfe1 <- lm(nlp_anger.cdc*100~-1+I(p_man*10)+log(n_100k.vac)+
                                 log_new_case_7ma+log_new_death_7ma+log_new_dose_7ma+
              I(Series_Complete_Pop_Pct/10)+I(Series_Complete_Pop_Pct_sd/10)+
                  as.factor(date)+as.factor(state),
                       data=data_panel)
vcov_nlp_anger.cdc_twfe1 <- vcovPL(lm_nlp_anger.cdc_twfe1, 
                                          cluster = data_panel$state,
                                          order.by = data_panel$date)

output_lm_nlp_anger.cdc_twfe1<-coeftest(lm_nlp_anger.cdc_twfe1, vcov = vcov_nlp_anger.cdc_twfe1)
output_lm_nlp_anger.cdc_twfe1<-output_lm_nlp_anger.cdc_twfe1[str_detect(row.names(output_lm_nlp_anger.cdc_twfe1),"^(?!as)"),]
output_lm_nlp_anger.cdc_twfe1
```

Figure 5
```{r}
data.frame(dv="Rule-Based\nNegativity",
           iv=factor(names(output_lm_vader_neg.cdc_twfe1[,1]),
                         levels = names(output_lm_vader_neg.cdc_twfe1[,1]),
                         labels = c("Mandate tweets%\n(Every 10%)",
                                    "log n\nvaccine tweets",
                                    "log new case",
                                    "log new death",
                                    "log new dose",
                                    "Vaccine coverage",
                                    "Vaccine coverage sd")),
           estimate=output_lm_vader_neg.cdc_twfe1[,1],
           lb95=output_lm_vader_neg.cdc_twfe1[,1]+qnorm(0.025)*output_lm_vader_neg.cdc_twfe1[,2],
           ub95=output_lm_vader_neg.cdc_twfe1[,1]+qnorm(0.975)*output_lm_vader_neg.cdc_twfe1[,2],
           lb90=output_lm_vader_neg.cdc_twfe1[,1]+qnorm(0.05)*output_lm_vader_neg.cdc_twfe1[,2],
           ub90=output_lm_vader_neg.cdc_twfe1[,1]+qnorm(0.95)*output_lm_vader_neg.cdc_twfe1[,2])%>%
  bind_rows(data.frame(dv="ML-based\nNegativity",
           iv=factor(names(output_lm_nlp_neg.cdc_twfe1[,1]),
                         levels = names(output_lm_nlp_neg.cdc_twfe1[,1]),
                         labels = c("Mandate tweets%\n(Every 10%)",
                                    "log n\nvaccine tweets",
                                    "log new case",
                                    "log new death",
                                    "log new dose",
                                    "Vaccine coverage",
                                    "Vaccine coverage sd")),
           estimate=output_lm_nlp_neg.cdc_twfe1[,1],
           lb95=output_lm_nlp_neg.cdc_twfe1[,1]+qnorm(0.025)*output_lm_nlp_neg.cdc_twfe1[,2],
           ub95=output_lm_nlp_neg.cdc_twfe1[,1]+qnorm(0.975)*output_lm_nlp_neg.cdc_twfe1[,2],
           lb90=output_lm_nlp_neg.cdc_twfe1[,1]+qnorm(0.05)*output_lm_nlp_neg.cdc_twfe1[,2],
           ub90=output_lm_nlp_neg.cdc_twfe1[,1]+qnorm(0.95)*output_lm_nlp_neg.cdc_twfe1[,2]))%>%
  bind_rows(data.frame(dv="Dictionary-based\nAnger",
           iv=factor(names(output_lm_NRC_anger.cdc_twfe1[,1]),
                         levels = names(output_lm_NRC_anger.cdc_twfe1[,1]),
                         labels = c("Mandate tweets%\n(Every 10%)",
                                    "log n\nvaccine tweets",
                                    "log new case",
                                    "log new death",
                                    "log new dose",
                                    "Vaccine coverage",
                                    "Vaccine coverage sd")),
           estimate=output_lm_NRC_anger.cdc_twfe1[,1],
           lb95=output_lm_NRC_anger.cdc_twfe1[,1]+qnorm(0.025)*output_lm_NRC_anger.cdc_twfe1[,2],
           ub95=output_lm_NRC_anger.cdc_twfe1[,1]+qnorm(0.975)*output_lm_NRC_anger.cdc_twfe1[,2],
           lb90=output_lm_NRC_anger.cdc_twfe1[,1]+qnorm(0.05)*output_lm_NRC_anger.cdc_twfe1[,2],
           ub90=output_lm_NRC_anger.cdc_twfe1[,1]+qnorm(0.95)*output_lm_NRC_anger.cdc_twfe1[,2]))%>%
  bind_rows(data.frame(dv="ML-based\nAnger",
           iv=factor(names(output_lm_nlp_anger.cdc_twfe1[,1]),
                         levels = names(output_lm_nlp_anger.cdc_twfe1[,1]),
                         labels = c("Mandate tweets%\n(Every 10%)",
                                    "log n\nvaccine tweets",
                                    "log new case",
                                    "log new death",
                                    "log new dose",
                                    "Vaccine coverage",
                                    "Vaccine coverage sd")),
           estimate=output_lm_nlp_anger.cdc_twfe1[,1],
           lb95=output_lm_nlp_anger.cdc_twfe1[,1]+qnorm(0.025)*output_lm_nlp_anger.cdc_twfe1[,2],
           ub95=output_lm_nlp_anger.cdc_twfe1[,1]+qnorm(0.975)*output_lm_nlp_anger.cdc_twfe1[,2],
           lb90=output_lm_nlp_anger.cdc_twfe1[,1]+qnorm(0.05)*output_lm_nlp_anger.cdc_twfe1[,2],
           ub90=output_lm_nlp_anger.cdc_twfe1[,1]+qnorm(0.95)*output_lm_nlp_anger.cdc_twfe1[,2]))%>%
  mutate(dv=factor(dv,levels=c("Rule-Based\nNegativity","ML-based\nNegativity","Dictionary-based\nAnger","ML-based\nAnger")),
         iv=fct_rev(iv))%>%
  ggplot(aes(iv, estimate))+
  geom_point()+
  geom_text(aes(label=round(estimate,3)),nudge_x = 0.4,size=2.5)+
  geom_errorbar(aes(ymin = lb95, ymax = ub95), width=.2, linewidth=0.3)+
  geom_errorbar(aes(ymin = lb90, ymax = ub90), width=.4, linewidth=0.3)+
  geom_hline(yintercept = 0,
             colour = gray(1/2), lty = 2)+
  scale_x_discrete(name="Independent Variables")+
  facet_wrap(~dv,scales="free")+
  scale_y_continuous(name="Effects")+
  coord_flip()+
  theme_bw(base_size=9)
ggsave("p_panel_cdc.tiff",width=1920,height=1280,units = "px")
```